#ifndef NAVIERSTOKESSOLVER_H
#define NAVIERSTOKESSOLVER_H

#include "stdafx.h"
#include "SceneObj.h"
#include "Utility.h"
#include "GridField.h"
#include "MarchingHeap.h"

class TwowayCoupler;
class DecoScene;
class CSRMatrix;
class FLuidSparseMatrix;

#define PARTICLE_PER_COORD 2
#define PARTICLE_PER_CELL (PARTICLE_PER_COORD * PARTICLE_PER_COORD * PARTICLE_PER_COORD)

#define RHO_WATER 1000.0
#define VISCOSITY_WATER 0.000691
static DOUBLE unifiedViscosity = VISCOSITY_WATER / RHO_WATER;
static const INT blockXWidth = 16;
static const INT blockYWidth = 16;
static const INT blockZWidth = 16;

enum FaceDir
{
	FD_UAxis,
	FD_VAxis,
	FD_WAxis
};


class NSSolver
{
private:
	DOUBLE m_viscocity;
	DOUBLE m_tension;
	vector3 m_force;
	INT m_blockSizeX; //num blocks in X direction
	INT m_blockSizeY;
	INT m_blockSizeZ;
	INT m_showLayer;
	DecoScene* m_scene;
	//list<Particle> m_markerParticles;
	list<Index3> m_blockToUpdate;
	GridField* m_grid;
	MarchingHeap heap;
	GridField resultField;
	GridField errorField;
	GridField resultField2;
	GridField tmpMinusVelField;
	DOUBLE*** divergenceVelocity;
	DOUBLE*** ADiag;
	DOUBLE*** APlusi;
	DOUBLE*** APlusj;
	DOUBLE*** APlusk;
	BOOL*** tmpDomain;
	//DOUBLE*** p;
	DOUBLE*** s;
	DOUBLE*** z;
	DOUBLE*** precon;
	DOUBLE*** r;
	DOUBLE*** q;
	DOUBLE*** tmpU;
	DOUBLE*** tmpV;
	DOUBLE*** tmpW;
	TwowayCoupler* mCoupler;
	bool mbMatrixFormed;

public:
	BYTE*** solidElements;

private:
	void destroyAllArrays();
	void allocateAllArrays(INT numX, INT numY, INT numZ);
	void addForce();
	void addForce(GridField* grid);
	vector3 semiLagrange(const GridField& valueGrid, const vector3& pos, const vector3& velocity, DOUBLE timeStep);
	void convectWithoutBFECC(const GridField& velocityField, const GridField& valueField, GridField& resultField, const vector<Index3>& indicesU, const vector<Index3>& indicesV, const vector<Index3>& indicesW);
	void convection(GridField* grid);
	void convection();
	void projection(GridField* grid, DOUBLE*** divergenceVelocity, DOUBLE*** ADiag, DOUBLE*** APlusi, DOUBLE*** APlusj, DOUBLE*** APlusk, BOOL*** tmpDomain);
	void projection();
	void extrapolateVelocity();
	void calculateVorticity();
	void calculateVorticity(GridField* grid);
	DOUBLE getPressure(GridField* grid, UINT i, UINT j, UINT k);
	DOUBLE getU(GridField* grid, UINT i, UINT j, UINT k);
	DOUBLE getV(GridField* grid, UINT i, UINT j, UINT k);
	DOUBLE getW(GridField* grid, UINT i, UINT j, UINT k);


	void initializeLevelSet(GridField* marchingGrid);
	DOUBLE fastSweepThroughOneElement(GridField* marchingGrid, INT i, INT j, INT k);
	void fastSweepAllDirs();
	void fastMarching(GridField* marchingGrid, DOUBLE marchingDist);
	void enforceSolidVeloctiy();
	void enforceSolidVeloctiy(GridField* grid);
	void calculateCloseElement(GridField* marchingGrid, LevelSetWithIndex& elem);

	void updateGridType(GridField* grid);
	void updateGridType();

	void formMatrixAndDivergence(GridField* grid, DOUBLE ***ADiag, DOUBLE ***APlusi, DOUBLE ***APlusj, DOUBLE ***APlusk, DOUBLE ***divergenceVelocity, vector<Index3>& indices, BOOL*** domain);
	BOOL residualApproximateZero(DOUBLE ***residual, const vector<Index3>& indices, INT numX, INT numY, INT numZ);
	void computePreconditioner(GridField* grid, DOUBLE ***ADiag, DOUBLE ***APlusi, DOUBLE ***APlusj, DOUBLE ***APlusk, DOUBLE ***precon, const vector<Index3>& indices);
	void applyPreconditioner(GridField* grid, DOUBLE ***ADiag, DOUBLE ***APlusi, DOUBLE ***APlusj, DOUBLE ***APlusk, DOUBLE ***precon, DOUBLE*** vec, DOUBLE*** q, const vector<Index3>& indices, DOUBLE ***result);
	void applyMatrixA(GridField* grid, DOUBLE ***ADiag, DOUBLE ***APlusi, DOUBLE ***APlusj, DOUBLE ***APlusk, DOUBLE*** vec, const vector<Index3>& indices, DOUBLE ***result);
	DOUBLE dotProduct(DOUBLE*** v1, DOUBLE*** v2, const vector<Index3>& indices);

	void solveForPressure(GridField* grid, DOUBLE ***ADiag, DOUBLE ***APlusi, DOUBLE ***APlusj, DOUBLE ***APlusk, DOUBLE ***divergenceVelocity, const vector<Index3>& indices, BOOL*** domain);
	void projectVelocity(GridField* grid, const vector<Index3>& indices, BOOL*** domain);
public:
	void checkIncompressible(GridField* grid);
	void FormPoissonMatrix();
	CSRMatrix* GetCSRMatrix();
	FLuidSparseMatrix* GetFluidMatrix();
	void SetGridField(GridField* gridField);
	void SetScene(DecoScene* scene);
	void UpdateFluidGridType(bool bRefillParticles = false);
	void Solve ();
	void ToggleShowGridLayer();
	void Render (DecoRenderInterface* RI, DecoLight** Lights, UINT numEffectiveLights);
	vector3 GetVelocity(const vector3& pos, BOOL fromCoarse = FALSE);	
	//void fillWithFluid();
	void fillAllGridWithFluid(INT airMargin, bool bRefillParticles);
	void fillAllGridWithFluid(GridField* grid, INT airMargin);
	void Serialize(DecoArchive& Ar);
	void Deserialize(DecoArchive& Ar);
	NSSolver();
	~NSSolver()
	{
		//m_markerParticles.clear();
		destroyAllArrays();
	}
	//static NSSolver* GetSingleton();
	//static void DestroySingleton();
	friend class NSSolverTet;
};


#endif
